转录组 ( Transcriptome ) 是指在特定环境或生理条件下的一个细胞或一个细胞群中所有 RNA 分子的集合,包括信使 RNA ( mRNA )、核糖体 RNA ( rRNA )、转运 RNA ( tRNA )及非编码 RNA ( ncRNA )。它有时被用来指所有的 RNA ,或仅指 mRNA ,指何种 RNA 取决于我们设计的特定实验,以下转录组分析仅针对 mRNA 进行研究。转录组不同于正常基因组的稳定性,它会随外界环境的变化而发生动态变化。因此,转录组可以反映不同实验条件和不同时间里基因表达的变化,但是无法获取转录衰减等 mRNA 降解现象的信息。根据实验设计获得不同的实验样本后,提取 mRNA 并通过二代测序获取 RNA 序列信息用于后续信息分析,本报告分析内容主要以转录组的差异基因分析为重点,并进行功能注释来进一步探究不同的实验设计或环境是如何影响样本的基因表达和功能的发挥
转录组测序一般需要经过 RNA 样品提取,RNA 检测,测序文库构建,上机测序等一系列实验环节,最终得到用于生物信息分析的测序数据。其中每一个环节都会影响测序数据的产出和质量,从而影响后续生物信息分析结果的准确性。镁伽科技重视每一个实验环节的准确性和可靠性,以确保得到高质量的测序数据。具体建库测序流程如下:
建库流程
高质量的RNA是整个项目成功的基础,为保证测序数据准确性,我们对 RNA 样品的检测主要包括以下几点:
样品检测合格之后,进行样品的文库构建,具体mRNA 捕获建库过程如下:
文库构建成功之后对文库进行检测,库检主要包括以下两方面:
库检合格后,把不同文库按照有效浓度及目标下机数据量的需求 pooling 后进行 Illumina PE150 测序。
目前 mRNA 转录组测序生物信息分析过程主要包括数据质控,参考基因组比对,基因表达量分析,差异基因分析,功能富集分析、差异基因的蛋白互作分析,变异检测、差异可变剪接分析、GSEA分析。分析流程如下图所示:
信息分析流程
测序获得的原始数据 ( Raw Data ) 存在一定比例低质量数据,会对后续信息分析造成很大干扰。高质量的测序数据是后续信息分析结果可靠性的前提,因此需要对原始数据进行严格的数据预处理和评估。我们使用自主研发的开源软件 fastp 1对原始数据进行处理,对低质量序列进行过滤,过滤内容包括如下部分:
测序过程本身存在机器错误的可能性,测序错误率分布检查可以反映测序数据的质量,序列信息中每个碱基的测序质量值保存在 fastq 文件中。如果测序错误率用 e 表示,Illumina 的碱基质量值用 Qphred 表示,则有:Qphred=-10log10(e) 。Illumina Casava 1.8 版本碱基识别与 Phred 分值之间的简明对应关系见下表。
| Phred 分值 | 不正确的碱基识别 | 碱基正确识别率 | Q.sorce |
|---|---|---|---|
| 10 | 1/10 | 90% | Q10 |
| 20 | 1/100 | 99% | Q20 |
| 30 | 1/1000 | 99.90% | Q30 |
| 40 | 1/10000 | 99.99% | Q40 |
选取STAR ( Spliced Transcripts Alignment to a Reference )软件2 将转录组测序 Reads 比对到参考基因组。STAR具有速度快、灵敏度和准确度高、内存消耗低等优势。Clean Reads 比对到参考基因组上获得比对特征信息并以 SAM 格式( Sequence Alignment/Map format ) 存储。SAM 序列比对格式标准由sanger制定 3 ,主要用于测序序列比对到基因组上信息特征的展示,BAM 文件是 SAM 文件的二进制格式,存储资源更小,常用于后续计算。基于比对结果使用 Samtools 对 BAM 文件进行排序、过滤重复 Reads 等处理,最后对序列在参考基因组上的比对结果进行统计和绘图展示。
数据下机经过质控处理后,将数据比对到参考基因组上,具体的比对结果统计见文件: results/qc/multiqc_report.html
表达量分析是差异分析的前提,表达量分析的实质是根据各样本的比对结果统计不同基因的 Reads Count 值,即比对到该基因的 Reads 数有多少,也是后续差异基因分析的基础信息,如果基因的表达量高,则其 Reads 覆盖数也会高。
| gene | N12 | N13 | N15 | T12_D | T12_NC | T13_D | T13_NC | T15_D | T15_NC |
|---|---|---|---|---|---|---|---|---|---|
| 5S_rRNA | 6 | 0 | 1 | 0 | 1 | 1 | 5 | 1 | 0 |
| 5_8S_rRNA | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 7SK | 1 | 3 | 2 | 4 | 0 | 3 | 0 | 0 | 0 |
| A1BG | 0 | 1 | 3 | 0 | 0 | 3 | 0 | 1 | 0 |
| A1BG-AS1 | 4 | 2 | 4 | 2 | 1 | 0 | 0 | 2 | 1 |
| A1CF | 1 | 2 | 1 | 1 | 5 | 0 | 0 | 1 | 3 |
TPM箱线图
样本TPM表达矩阵相关性热图
pca 散点图
差异分析 ( Difference analysis ) 是在一组样本中找出产生差异的原因及其对造成差异的影响程度。根据基因的定量分析信息可进行差异基因分析,筛选不同的实验条件下的功能基因。基因差异表达分析的输入数据为基因表达水平分析中得到的Reads Count数据,分析主要分为三部分:
| 类型 | 软件 | 标准化方法 | p值矫正方法 |
|---|---|---|---|
| 有生物学重复 | Deseq24 | 软件自设 | BH5 |
| 无生物学重复 | EdgeR6 | 软件自设 | BH |
基因表达具有时间和空间特异性,在两个不同条件下,表达水平存在显著差异,各比较组合间分析得到的差异表达基因统计如下表所示:
| gene | baseMean | log2FoldChange | lfcSE | pvalue | padj |
|---|---|---|---|---|---|
| LINC03061 | 1456.5710438 | -2.2579168 | 0.2199267 | 0.0000000 | 0 |
| FSTL1 | 320.3733667 | -3.2547143 | 0.3367156 | 0.0000000 | 0 |
| SCEL | 746.2850365 | -1.8190108 | 0.2019606 | 0.0000000 | 0 |
| REG1A | 495.0132112 | 5.2736542 | 0.6889077 | 0.0000000 | 0 |
| TMEM238L | 294.8064412 | 1.7754076 | 0.2359544 | 0.0000000 | 0 |
| SERPINB2 | 186.6521241 | -4.4830945 | 0.5878698 | 0.0000000 | 0 |
| DDX11L2 | 0.3177327 | 0.0000049 | 0.0089362 | 0.8024778 | 1 |
| WASH7P | 128.3431754 | 0.0001754 | 0.0089322 | 0.5105055 | 1 |
| MIR6859-1 | 2.5199039 | -0.0000302 | 0.0089360 | 0.6187630 | 1 |
| MIR1302-2HG | 0.2261730 | 0.0000000 | 0.0089362 | 1.0000000 | 1 |
| ENSG00000238009 | 0.3863901 | -0.0000056 | 0.0089362 | 0.7771276 | 1 |
| CICP27 | 3.0549363 | -0.0000006 | 0.0089360 | 0.9924111 | 1 |
差异基因表达火山图
基因表达数目韦恩图
通过富集分析,可以找到不同条件下的差异基因与哪些生物学功能或通路显著性相关,从而揭示和理解生物学过程的基本分子机制。 我们采用 clusterProfiler7 软件对差异基因集进行GO功能富集分析,KEGG通路富集分析等。富集分析基于超几何分布原理,其中差异基因集为差异显著分析所得差异基因并注释到GO或KEGG数据库的基因集,背景基因集为所有参与差异分析并注释到GO或KEGG数据库的基因集。富集分析结果是对每个差异比较组合的所有差异基因集进行富集。
富集分析原理图
GO(Gene Ontology)是描述基因功能的综合性数据库(http://www.geneontology.org/),可分为生物过程(Biological Process)和细胞组成(Cellular Component)分子功能(Molecular Function)三个部分。GO 功能显著性富集分析给出与基因集背景相比,在差异表达基因中显著富集的 GO 功能 条目,从而给出差异表达基因与哪些生物学功能显著相关。该分析首先把所有差异表达基因向 Gene Ontology 数据库的各个 term 映射,计算每个 term 的基因数目,然后找出与整个基因集背景相比,在差异表达基因中显著富集。GO富集分析方法为clusterProfiler的enricher函数,该方法使用超几何分布检验的方法获得显著富集的GO Term。统计被显著富集的各个GOterm中的基因数,GO富集以padj小于0.05为显著富集。
GO_BP
GO_MF
GO_CC
在生物体内,不同基因相互协调行使其生物学功能,通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。 KEGG(Kyoto Encyclopedia of Genes and Genomes)是有关Pathway的主要公共数据库,其中整合了基因组化学和系统功能信息等众多内容。信息分析时,我们将KEGG数据库按照动物,植物,真菌等进行分类,依据研究物种的种属选择相应的类别进行分析。Pathway显著性富集分析以KEGG Pathway为单位,应用超几何检验,找出差异基因相对于所有有注释的基因显著富集的pathway,KEGG通路富集同样以padj小于0.05作为显著性富集的阈值。
蛋白互作网络 ( Protein-Protein Interaction Networks ,PPI ) 是蛋白通过彼此之间的相互作用构成,参与基因表达调节、生物信号传递、细胞周期调控及能量物质代谢等生命过程的各个环节。我们使用 STRING ( https://string-db.org/ ) 蛋白质互作数据库中的互作关系进行差异基因的 PPI 分析,对于差异基因的蛋白互作网路数据文件可使用 StringDB 进行可视化展示8。
详细的差异基因蛋白互作网络见文件result/enrichment/ * / * /STRING_info.txt
可变剪接 ( Alternative Splicing ) 是指一个 mRNA 前体通过不同的剪接方式 ( 选择不同的剪接位点 ) 产生不同的 mRNA 剪接异构体,是真核生物基因和蛋白质数量有较大差异的重要原因。可变剪接可分为
可变剪接分析原理
每个比较组合得到的差异可变剪接统计如下表:
| EventType | EventTypeDescription | TotalEventsJC | TotalEventsJCEC | SignificantEventsJC | SigEventsJCSample1HigherInclusion | SigEventsJCSample2HigherInclusion | SignificantEventsJCEC | SigEventsJCECSample1HigherInclusion | SigEventsJCECSample2HigherInclusion |
|---|---|---|---|---|---|---|---|---|---|
| SE | skipped exon | 83845 | 83845 | 0 | 0 | 0 | 0 | 0 | 0 |
| A5SS | alternative 5’ splice sites | 9546 | 9546 | 0 | 0 | 0 | 0 | 0 | 0 |
| A3SS | alternative 3’ splice sites | 14512 | 14512 | 0 | 0 | 0 | 0 | 0 | 0 |
| MXE | mutually exclusive exons | 8716 | 8716 | 0 | 0 | 0 | 0 | 0 | 0 |
| RI | retained intron | 9881 | 9881 | 0 | 0 | 0 | 0 | 0 | 0 |
| SE | skipped exon | 85040 | 85040 | 0 | 0 | 0 | 0 | 0 | 0 |
| A5SS | alternative 5’ splice sites | 9613 | 9613 | 0 | 0 | 0 | 0 | 0 | 0 |
| A3SS | alternative 3’ splice sites | 14564 | 14564 | 0 | 0 | 0 | 0 | 0 | 0 |
| MXE | mutually exclusive exons | 8975 | 8975 | 0 | 0 | 0 | 0 | 0 | 0 |
| RI | retained intron | 9895 | 9895 | 0 | 0 | 0 | 0 | 0 | 0 |
GSEA分析(Gene Set Enrichment Analysis),基因集富集分析9,是目前非常常用的RNA表达分析手段。GSEA不需要指定明确的差异基因阈值,算法会根据实际数据的整体趋势,为我们提供一种合理的解决方法,即使在没有先验经验的情况下也能在表达谱整体层次上对多个基因进行分析,从而从数理统计层面把表达谱数据与生物学意义很好地衔接起来。目前分析物种只支持人、小鼠,最低条件满足3V3组合才能跑gsea分析。
下图展示了基因表达的可变性与细胞中平均表达水平的关系,红线显示的了基于此关系的期望值, 我们将用于后续的细胞亚群识别或沿轨迹排序细胞标记为黑点,而其他基于显示为灰点。
| X | ID | Description | setSize | enrichmentScore | NES | pvalue | p.adjust | qvalue | rank | leading_edge | core_enrichment |
|---|---|---|---|---|---|---|---|---|---|---|---|
| hsa04390 | hsa04390 | Hippo signaling pathway | 129 | -0.6189196 | -2.311827 | 0.0000004 | 0.0000234 | 0.0000107 | 1187 | tags=33%, list=8%, signal=31% | 7529/6788/7533/659/5518/10971/3398/23291/324/4771/4087/126374/7004/5515/7046/7532/2932/5501/896/5516/25937/26524/658/55233/7534/999/60/332/7477/8323/154796/84962/71/2736/1490/4092/1739/81029/23286/7476/595/894/3689 |
| hsa03010 | hsa03010 | Ribosome | 129 | -0.6128998 | -2.289342 | 0.0000008 | 0.0000234 | 0.0000107 | 3903 | tags=70%, list=26%, signal=52% | 63931/29074/6231/11224/9349/6147/6135/9045/6139/6132/6222/6156/6154/6169/219927/6157/6168/6165/29093/6142/7311/6188/6209/6183/6208/51121/6217/51264/6143/6161/6227/6155/6160/6166/51263/6228/6130/6230/6171/11222/6133/6152/6146/6234/6235/6137/6176/6129/6189/6218/55052/6122/51187/6136/6187/51318/6125/6210/6181/6150/6204/6134/51021/6138/6124/6201/6158/28998/6175/3921/64928/6232/6223/6164/6203/4736/29088/55173/6224/9553/23521/9801/200916/6193/64965/6144/6207/388524/6233/6206 |
| hsa04110 | hsa04110 | Cell cycle | 154 | -0.4910432 | -1.864148 | 0.0005038 | 0.0104127 | 0.0047732 | 4570 | tags=57%, list=30%, signal=40% | 9184/5519/10926/991/10735/23595/151648/4176/27085/1387/23047/113130/4089/51434/25/81620/3065/993/4998/4171/8317/9978/9088/890/8318/9134/1019/8697/4174/57082/8454/4175/4173/51529/1031/701/2810/5347/1871/983/891/9133/10393/5885/64682/7531/1870/6502/3066/7029/4085/1869/5001/5111/8379/6500/5591/23244/9319/10403/5925/5527/7529/545/7533/23063/5518/990/10971/1111/5529/157570/4087/5515/7532/25847/5933/2932/896/5516/7027/7534/7157/10274/1021/1029/595/894 |
| hsa05132 | hsa05132 | Salmonella infection | 213 | -0.4463716 | -1.772945 | 0.0006940 | 0.0107566 | 0.0049309 | 3402 | tags=40%, list=22%, signal=32% | 6416/9367/1499/1432/10376/1781/5291/26999/6188/4074/1434/203068/4790/5605/1639/5868/302/5216/8767/10552/5585/5606/51143/5878/10392/140735/207/5287/836/10121/65082/6885/3840/10787/84516/3265/3836/5898/375/10093/5601/10006/6500/81873/382/6281/27072/55207/1778/84790/3799/55770/5594/3839/8976/10096/998/387/8772/1147/5420/11258/10671/9842/10097/10640/112574/54205/4615/10627/8655/103910/55860/840/3654/60/55845/79026/5288/10092/71/9265/7846/3606/2316/10398 |
| hsa04360 | hsa04360 | Axon guidance | 144 | -0.4904220 | -1.846516 | 0.0012931 | 0.0111653 | 0.0051182 | 1395 | tags=24%, list=9%, signal=22% | 3688/5534/5594/818/4233/23365/659/8440/5361/998/387/5295/1072/2773/23654/3845/5062/10627/1946/2932/1948/658/103910/4775/29984/4893/5781/1969/1808/10371/81029/8503/56963/10512/10398 |
| hsa04530 | hsa04530 | Tight junction | 142 | -0.4852686 | -1.826337 | 0.0014407 | 0.0111653 | 0.0051182 | 1689 | tags=30%, list=11%, signal=27% | 5111/10093/5601/81873/7074/3996/3688/84790/1366/10096/998/5518/387/153562/81/11346/10097/4771/5515/10627/5516/1364/100506658/7082/103910/60/154810/137075/79784/51762/10092/154796/4627/71/9076/7846/7430/1739/51421/5581/595/10398 |
| hsa04510 | hsa04510 | Focal adhesion | 186 | -0.5988323 | -2.061757 | 0.0000000 | 0.0000000 | 0.0000000 | 3027 | tags=31%, list=12%, signal=27% | 3918/5595/4659/3655/5291/8394/1398/71/7408/29780/10627/6655/6093/3611/3694/10398/5747/5500/3673/3678/2013/3696/5604/4233/5501/5170/8395/673/5829/56034/2909/3691/5599/5881/3675/2317/3909/3913/5906/998/3688/7414/6654/329/2932/331/4660/595/824/2316/1956/5601/2064/340156/5908/394/7791 |
| hsa04810 | hsa04810 | Regulation of actin cytoskeleton | 206 | -0.5387199 | -1.869518 | 0.0000000 | 0.0000006 | 0.0000003 | 3757 | tags=34%, list=15%, signal=29% | 4628/1793/5216/26999/27006/9459/10000/8826/2254/10097/5962/10109/5595/4659/3655/5291/85464/8394/1398/71/10152/10096/10627/6655/6093/3694/54434/3984/10092/10398/3845/5747/10163/5500/3673/3678/2934/3696/5604/8936/55740/6548/28964/5501/8395/673/5829/56034/1730/2909/3691/81624/5881/3675/324/998/3688/7414/6654/10297/4660/55970/4629/1956/23365/55845/200576/340156/79784/4627/85477 |
| hsa04820 | hsa04820 | Cytoskeleton in muscle cells | 203 | -0.5145934 | -1.786836 | 0.0000008 | 0.0000230 | 0.0000118 | 4996 | tags=43%, list=20%, signal=35% | 25802/58/482/2026/57731/3685/2274/6443/287/7139/1158/6385/7170/29765/6382/23224/3674/91977/89/83660/3908/3695/93661/2010/64236/7273/9499/9672/481/1829/23002/1292/3672/6712/7111/1277/57644/2192/4628/832/1605/29767/9260/3339/6641/1824/8557/8910/79188/3728/3655/6709/22989/478/5318/9124/6711/71/7168/476/271/161176/3694/10611/10398/3673/3678/3696/1466/1730/829/3691/81624/3675/3688/7414/7134/288/4629/6645/29766/2275/7171/1832/7791/79784/4627 |
| hsa04530 | hsa04530 | Tight junction | 158 | -0.5167376 | -1.742798 | 0.0000097 | 0.0002095 | 0.0001077 | 3782 | tags=39%, list=15%, signal=34% | 5566/4628/4217/57826/10097/51422/5962/5520/9080/10109/55114/4771/9076/861/9414/2017/9693/51735/5581/51421/71/3993/7408/10207/10096/9368/10627/6093/10092/10398/54566/9223/50848/5562/11346/9181/5563/55844/100506658/84952/5584/57530/56288/5599/92359/123720/5906/154796/998/3688/4301/137075/595/1739/4629/5601/2064/4218/79784/5590/4627/51762 |
| hsa05165 | hsa05165 | Human papillomavirus infection | 301 | -0.4410204 | -1.587632 | 0.0000145 | 0.0002497 | 0.0001284 | 3512 | tags=33%, list=14%, signal=29% | 4855/8324/3716/1277/6772/245973/1027/3516/9134/5566/4854/8325/10000/6199/3455/5734/83439/1387/84441/1385/55012/8717/894/6009/3280/148022/5520/3918/10379/1857/5315/5595/182/4600/7098/3655/64764/5291/7976/472/5528/535/6773/528/5663/3993/54626/10207/6198/1017/9368/1021/6655/51382/5527/3694/6934/3845/5747/9223/3673/3678/3696/5604/5610/23352/55844/3066/7132/5829/9794/29110/5584/3691/56288/7476/5529/92359/5525/3675/3909/2033/3913/324/998/3688/8638/3065/8772/6654/10297/2932/595/3454/1739/1956/5523/7855/5590/55534 |
| hsa04144 | hsa04144 | Endocytosis | 243 | -0.4614140 | -1.630226 | 0.0000177 | 0.0002538 | 0.0001305 | 3471 | tags=33%, list=14%, signal=28% | 1211/1785/11031/10097/3303/9897/22841/60682/10109/3482/58513/9525/2060/50807/55040/23624/4087/116985/10564/378/10193/9765/8411/7048/8394/8724/84440/10096/137492/11267/387680/55048/10015/10092/57403/10617/382/7879/3304/10565/28964/11059/55616/51534/116987/8395/23325/57132/30844/5584/79720/829/8853/56288/9372/9135/408/9815/3949/6642/55680/998/80223/26286/80230/64744/83737/9265/5868/9267/116984/1956/30846/4218/64145/5869/5590/8766/9101 |
组织是由不同谱系和亚型的细胞类型组成的复杂环境,每种细胞都有自己独特的转录组。因此,批量转录组分析是细胞类型特异性基因表达的总和加权的细胞类型比例在给定的样本。去卷积的基因表达谱允许重建组织的细胞组成。xCell 是一个强大的计算方法,转换基因表达谱为丰富分数的免疫和基质细胞类型跨样本。
免疫浸润分数热图
| 分析内容 | 软件 | 版本 |
|---|---|---|
| 差异分析 | deseq2 | 1.38 |
| 差异分析 | edgeR | 3.40.2 |
| 富集分析 | clusterProfiler | 4.10.0 |
| 剪接分析 | rmats | 4.3.0 |
| 数据质控 | rseqc | 5.0 |
| 免疫浸润分析 | xcell | 1.3 |
Shifu, C. , Yanqing, Z. , Yaru, C. , & Jia, G. (2018). Fastp: an ultra-fast all-in-one fastq preprocessor. Bioinformatics, 34(17), i884-i890.↩︎
Dobin, A., et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics (Oxford, England), 29(1), 15–21. https://doi.org/10.1093/bioinformatics/bts635↩︎
Li, H. E. A. (2009). The sequence alignment / map ( sam ) format. Bioinformatics, 25(1 Pt 2), 1653-4.↩︎
Love, M. I. , Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol, 15(12), 550.↩︎
Benjamini, Y. , Hochberg, Y. (1995). Controlling the False Discovery Rate - a Practical and Powerful Approach to Multiple Testing. J Roy Stat Soc B, 57(1), 289-300.↩︎
Robinson, MD. , McCarthy, DJ. , Smyth, GK. (2010). edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics, 26, 139–140.↩︎
Yu, G. , Wang, L. G. , Han, Y. , & He, Q. Y. (2012). Clusterprofiler: an r package for comparing biological themes among gene clusters. Omics-a Journal of Integrative Biology, 16(5), 284-287.↩︎
Szklarczyk D, Gable AL, Nastou KC, Lyon D, Kirsch R, Pyysalo S, Doncheva NT, Legeay M, Fang T, Bork P, Jensen LJ, von Mering C (2021). “The STRING database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets.” Nucleic Acids Research (Database issue), 49.↩︎
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102:15545–15550.↩︎